/*
 * Matrix.cc
 *
 *  Created on: Sep 22, 2012
 *      Author: HUNG
 */

#include "Matrix.h"
#include <string.h>
#include <stdlib.h>
#include <omnetpp.h>

Matrix::Matrix()
{
    this->row = 0;
    this->col = 0;
    this->data = NULL;
}
Matrix::~Matrix()
{
    this->row = 0;
    this->col = 0;
    this->data = NULL;
}

Matrix::Matrix(int _row, int _col)
{
    this->row = _row;
    this->col = _col;
    this->data = InitArray(this->row, this->col);
}

Matrix::Matrix(int _row, int _col, double **A)
{
    this->row = _row;
    this->col = _col;
    this->data = InitArray(this->row, this->col);
    for(int i = 0; i < this->row; i++)
        for(int j = 0; j < this->col; j++)
            this->data[i][j] = A[i][j];
}

Matrix::Matrix(int _row, int _col, int isZezoOrIndentity)
{
    this->row = _row;
    this->col = _col;
    this->data = InitArray(this->row, this->col);

    if(isZezoOrIndentity == 0) // is Zezo Matrix
    {
        for(int i = 0; i < this->row; i++)
            for(int j = 0; j < this->col; j++)
                this->data[i][j] = 0;
    }
    else if(isZezoOrIndentity == 1) //is Indentity Matrix
    {
        for(int i = 0; i < this->row; i++)
            for(int j = 0; j < this->col; j++)
            {
                if(i == j) this->data[i][j] = 1;
                else this->data[i][j] = 0;
            }
    }
    else return;
}

Matrix MatrixZezo(int _row, int _col)
{
    Matrix MT = Matrix(_row, _col);
    for(int i = 0; i < MT.row; i++)
        for(int j = 0; j < MT.col; j++)
            MT.data[i][j] = 0;
    return MT;
}

Matrix MatrixIndentity(int _row, int _col)
{
    Matrix MT = Matrix(_row, _col);
    for(int i = 0; i < MT.row; i++)
        for(int j = 0; j < MT.col; j++)
        {
            if(i == j) MT.data[i][j] = 1;
            else MT.data[i][j] = 0;
        };
    return MT;
}
int Matrix::getRow()
{
    return this->row;
}

int Matrix::getCol()
{
    return this->col;
}

void Matrix::setData(int i, int j, double _data)
{
    this->data[i][j] = _data;
}

double Matrix::getData(int i, int j)
{
    return this->data[i][j];
}

Matrix Matrix::operator+(const Matrix &MT)
{
    if(this->row != MT.row || this->col != MT.col) return Matrix();
    else
    {
        Matrix X = Matrix(this->row, this->col);
        for(int r = 0; r < this->row; r++)
            for(int c = 0; c < this->col; c++)
                X.data[r][c] = this->data[r][c] + MT.data[r][c];
        return X;
    }
}

Matrix Matrix::operator-(const Matrix &MT)
{
    if(this->row != MT.row || this->col != MT.col) return Matrix();
    else
    {
        Matrix X = Matrix(this->row, this->col);
        for(int r = 0; r < this->row; r++)
            for(int c = 0; c < this->col; c++)
                X.data[r][c] = this->data[r][c] - MT.data[r][c];
        return X;
    }
}


Matrix Matrix::operator*(const Matrix &MT)
{
    if(this->col != MT.row) return Matrix();
    else
    {
        Matrix X = Matrix(this->row, MT.col);

        for(int r = 0; r < X.row; r++)
            for(int c = 0; c < X.col; c++)
            {
                double temp = 0;
                for(int k = 0; k < this->col; k++)
                    temp += this->data[r][k] * MT.data[k][c];
                X.data[r][c] = temp;
            };
        return X;
    }
}

Matrix Matrix::operator*(const double &a)
{
    Matrix MT = Matrix(this->row, this->col, this->data);
    for(int r = 0; r < MT.row; r++)
        for(int c = 0; c < MT.col; c++)
        {
            MT.data[r][c] *= a;
        };
    return MT;
}
Matrix Matrix::transpose()
{
    Matrix X(this->col, this->row);
    for(int r = 0; r < X.row; r++)
        for(int c = 0; c < X.col; c++)
            X.data[r][c] = this->data[c][r];
    return X;
}

double det(Matrix &MT)
{
    if(MT.row != MT.col)
    {
        exit(0);
    }
    int N = MT.row;
    double **A = InitArray(N, N);

    for(int i = 0; i < N; i++)
        for(int j = 0; j < N; j++)
            A[i][j] = MT.data[i][j];

    double Det = 1;
    int k, sign = 1;
    for (int i = 0; i < N - 1; i++)
    {
        if (A[i][i] == 0) // Nếu gặp phần tử trên đường chéo chính bằng 0 thì tìm hàng khác để đổi
        {
            k = i+1;
            while (k < N && A[k][i] == 0) k++;
            if (k == N) return 0;
            for (int j = i; j < N; j++)// Đổi hàng i với hàng k
            {
                swap(A[i][j], A[k][j]);
            }
            sign = -sign; // Định thức đổi dấu đó
        }
        for (int j = i + 1; j < N; j++)// Biển đổi để các phần tử cùng cột ở hàng dưới bằng 0 đó
        {
            double temp = -A[j][i]  / A[i][i];
            for (int k = i + 1; k < N; k++)
                A[j][k] += temp * A[i][k]; // Nhân hàng i với (-a[j][i]/a[i][i]) rồi cộng vào hàng j
        }
        Det *= A[i][i]; // Tính dần det(a)
    }

    Det *= A[N-1][N-1] * sign; // Nhân với phần tử cuối chéo chính nữa là xong
    return Det;
}

Matrix Matrix::inverse()
{
    Matrix MT = Matrix(this->row, this->col, this->data);
    if(det(MT) == 0) return Matrix();
    else
    {
        int N = MT.row;
        int k;
        double **A = InitArray(N, N);
        for(int r = 0; r < N; r++)
            for(int c = 0; c < N; c++)
                A[r][c] = MT.data[r][c];

        Matrix B = MatrixIndentity(N, N); //Ma tran don vi

        for(int i = 0; i < N; i++) //Xu ly tu hang dau den hang cuoi
        {
            if(A[i][i] == 0)
            {
                k = i + 1;
                while(k < N && A[k][i] == 0) k++;
                for(int j = 0; j < N; j++) //Doi hang i voi k cua ma tran A, B
                {
                    swap(A[i][j], A[k][j]);
                    swap(B.data[i][j], B.data[k][j]);
                }
            }

            double temp = A[i][i];
            for(int j = i; j < N; j++)
                A[i][j] /= temp;
            for(int j = 0; j < N; j++)
                B.data[i][j] /= temp;

            for(int j = i + 1; j < N; j++)
            {
                double temp = A[j][i];
                for(int k = i; k < N; k++)
                    A[j][k] -= temp * A[i][k];
                for(int k = 0; k < N; k++)
                    B.data[j][k] -= temp * B.data[i][k];
            }
        }//Ket thuc qua trinh khu Gauss

        for(int i = N - 1; i > 0; i--) //Bat dau qua trinh Jordan
            for(int j = i - 1; j >= 0; j--)
            {
                double temp = A[j][i];
                for(int k = N - 1; k >= i; k--)
                    A[j][k] -= temp * A[i][k];
                for(int k = 0; k < N; k++)
                    B.data[j][k] -= temp * B.data[i][k];
            };
        return B;
    }
}

double** InitArray(int _row, int _col)
{
    double** A = new double *[_row];
    for(int i = 0; i < _row; i++)
        A[i] = new double[_col];
    return A;
}

void swap(double &x, double &y)
{
    double temp = x;
    x = y;
    y = temp;
}

void Matrix::Output()
{
    if(row == 0 || col == 0)
    {
        ev <<"\n NULL ";
    }
    for(int r = 0; r < row; r++)
        for(int c = 0; c < col; c++)
        {
            if( c == 0) ev << "\n";
            ev << "    " << data[r][c];
        };
}
